#include <stdio.h>
#include <graphics.h>
#include <math.h>


double*solvematrix(double**M,double*P,int N)
{
    double*X=new double[N];
 
    int i,j,m; 
    double K,S=0;
    double A[N][N];  
    double B[N]; 
    
    for (j=0; j<N; j++) 
    {
      for (i=0; i<N; i++) 
      {
        A[i][j]=*(((double*)M)+i*200+j);
      }
      B[j]=P[j];
    }

    for (i=0; i<=N-2; i++)
    {
      for (j=i+1;j<=N-1; j++)
      {
        K=A[i][i]/A[i][j];
        for (m=i; m<=N-1; m++) 
        {
          A[m][j]=A[m][j]*K-A[m][i];
        }
        B[j]=B[j]*K-B[i];
      }
    }
    
    for (j=N-1; j>=0; j--) 
    {
      S=0;
      for (i=N-1; i>j; i--) 
      {
        S+=A[i][j]*(*(X+i));
      }
      *(X+j)=(B[j]-S)/A[j][j];
    }
    
    return X;
}

int main() 
{
    int N=4; 
            
    double A[200][200];
    double B[200];
    double*X=new double[N];

    double Norm;
    int Step=0;
    double e=0.01;
    int i,j=0;
    
    initwindow(400, 300, "Norm(N)");
    setcolor(BLACK); 
    setbkcolor(WHITE);

    setfillstyle(1,RED);
    cleardevice();
    
    int x0=100, y0=200;   
    int x1=300, y1=100;   

    moveto(x0,y0);
    lineto(x0,y1);
    moveto(x0,y0);
    lineto(x1,y0);

    moveto(x0-3,y1+5);
    lineto(x0,y1);
    moveto(x0+3,y1+5);
    lineto(x0,y1);
    
    moveto(x1,y0);
    lineto(x1-5,y0-3);
    moveto(x1-5,y0+3);
    lineto(x1,y0);

    X=solvematrix((double**)A, (double*)B, N);
    
    for (i=0; i<=N-1; i++)
    {
      printf("x%d = %f; ",i + 1,*(X+i));
    }
    printf("\n");
    
    N=1;
    while(true==true)
    {
      N++;
      double*X1=new double[N];
      printf("N = %d\n\n", N);
      for (j=0; j<=N-1; j++)
      {
        for (i=0; i<=N-1; i++)
        {
          if(i==j)  
          {
            A[i][j]=1+e;
          }
          else if(i<j)
          {
            A[i][j]=e;
          }
          else if(i>j)
          {
            A[i][j]=1-e;
          }
        }
        if(j==N-1)
        {
          B[j]=1+e;
        }
        else
        {
          B[j]=-1;
        }
      }   
      
      for (j=0; j<=N-1; j++) 
      {
        for (i=0; i<=N-1; i++) 
        {
          printf("%f ", A[i][j]);
        }
          printf("-- %f\n", B[j]);
      }   
      printf("\n");

      X1 = solvematrix((double**)A, (double*)B,N);
      for (i=0; i<=N-1; i++)
      {
        printf("x%d = %.15f\n",i+1,*(X1+i));
        Norm+=X1[i]*X1[i];
      }
      printf("\n");
        
      Norm=pow(Norm, 0.5)/N;
      printf("Norm = %.15f\n\n\n", Norm);
      fillellipse(101+Step*5,200-(int)(Norm*100),1,1);
      Step++;
      delete X1;
      if(Step>=16)
      {
        break;   
      }
    }  
    
    getchar();
    closegraph();
	return 0;
}
